/*                     N U R B _ T R I M . C
 * BRL-CAD
 *
 * Copyright (c) 1990-2025 United States Government as represented by
 * the U.S. Army Research Laboratory.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * version 2.1 as published by the Free Software Foundation.
 *
 * This library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this file; see the file named COPYING for more
 * information.
 */
/** @addtogroup */
/** @{ */
/** @file primitives/bspline/nurb_trim.c
 *
 * Trimming curve routines.
 *
 */
/** @} */

#include "common.h"

#include <math.h>
#include "bio.h"

#include "vmath.h"
#include "nmg.h"

extern void rt_clip_cnurb(struct bu_list *plist, struct edge_g_cnurb *crv, fastf_t u, fastf_t v);

struct _interior_line {
    int axis;
    fastf_t o_dist;
};


#define QUAD1 0
#define QUAD2 1
#define QUAD3 2
#define QUAD4 3


#define TRIM_OUT 0
#define TRIM_IN 1
#define TRIM_ON 2

/* The following defines need to be 0, 2, 3 in order to drive the quad
 * table and determine the appropriate case for processing the
 * trimming curve.
 */

#define CASE_A 0
#define CASE_B 2
#define CASE_C 3

static int quad_table[16]  = {
    /* A = 0, B = 2, C = 3 */
    0, 0, 0, 0, 0, 3, 0, 3, 0, 2, 3, 3, 0, 3, 3, 3
};


/**
 * This routine determines what quadrants the trimming curves lies in,
 * It then uses a table look up to determine the whether its CASE{A,
 * B, C}, One difference from the paper is the fact that if any of the
 * points lie on the axis of the u, v quadrant system then the axis is
 * only in either Quadrant 1 or Quadrant 2 and not q3 or q4. This
 * handles the case of endpoint problems correctly.
 */
int
rt_trim_case(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
{
    int quadrant;
    int qstats;
    fastf_t * pts;
    int coords, rat;
    int i;

    qstats = 0;

    coords = RT_NURB_EXTRACT_COORDS(trim->pt_type);
    pts = trim->ctl_points;
    rat = RT_NURB_IS_PT_RATIONAL(trim->pt_type);

    /* Handle rational specially since we need to divide the rational
     * portion.
     */

    if (rat) {
	for (i = 0; i < trim->c_size; i++) {
	    if (pts[0]/pts[2] > u)
		quadrant = (pts[1]/pts[2] >= v)? QUAD1:QUAD4;
	    else
		quadrant = (pts[1]/pts[2] >= v)? QUAD2:QUAD3;

	    qstats |= (1 << quadrant);
	    pts += coords;
	}
    } else {
	for (i = 0; i < trim->c_size; i++) {
	    if (pts[0] > u)
		quadrant = (pts[1] >= v)? QUAD1:QUAD4;
	    else
		quadrant = (pts[1] >= v)? QUAD2:QUAD3;

	    qstats |= (1 << quadrant);
	    pts += coords;
	}
    }

    return quad_table[qstats];	/* return the special case of the curve */
}


/**
 * Process Case B curves.
 *
 * If the two endpoints of the curve lie in different quadrants than
 * the axis crosses the curve an odd number of times
 * (TRIM_IN). Otherwise the curve crosses the u, v axis a even number
 * of times (TRIM_OUT).  No further processing is required.
 */
int
rt_process_caseb(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
{
    int q1, q2;
    fastf_t * pts;
    int rat;

    rat = RT_NURB_IS_PT_RATIONAL(trim->pt_type);

    pts = trim->ctl_points;

    if (rat) {
	if (pts[0]/pts[2] > u) q1 = (pts[1]/pts[2] >= v)?QUAD1:QUAD4;
	else q1 = (pts[1]/pts[2] >= v)?QUAD2:QUAD3;


	pts = trim->ctl_points + RT_NURB_EXTRACT_COORDS(trim->pt_type) *
	    (trim->c_size -1);
	if (pts[0]/pts[2] > u) q2 = (pts[1]/pts[2] >= v)?QUAD1:QUAD4;
	else q2 = (pts[1]/pts[2] >= v)?QUAD2:QUAD3;

    } else {
	if (pts[0] > u) q1 = (pts[1] >= v)?QUAD1:QUAD4;
	else q1 = (pts[1] >= v)?QUAD2:QUAD3;


	pts = trim->ctl_points +
	    RT_NURB_EXTRACT_COORDS(trim->pt_type) 	*
	    (trim->c_size -1);
	if (pts[0] > u) q2 = (pts[1] >= v)?QUAD1:QUAD4;
	else q2 = (pts[1] >= v)?QUAD2:QUAD3;
    }

    if (q1 != q2)
	return TRIM_IN;
    else
	return TRIM_OUT;

}


/**
 * Only check end points of the curve
 */
int
rt_nurb_uv_dist(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
{

    fastf_t dist;
    fastf_t * ptr;
    int coords;
    int rat;
    fastf_t u2, v2;

    ptr = trim->ctl_points;
    coords = RT_NURB_EXTRACT_COORDS(trim->pt_type);
    rat = RT_NURB_IS_PT_RATIONAL(trim->pt_type);

    if (rat) {
	u2 = ptr[0]/ptr[2] - u; u2 *= u2;
	v2 = ptr[1]/ptr[2] - v; v2 *= v2;
    } else {
	u2 = ptr[0] - u; u2 *= u2;
	v2 = ptr[1] - v; v2 *= v2;
    }

    dist = sqrt(u2 + v2);
    if (NEAR_ZERO(dist, 1.0e-4))
	return TRIM_ON;

    ptr = trim->ctl_points + coords * (trim->c_size -1);

    if (rat) {
	u2 = ptr[0]/ptr[2] - u; u2 *= u2;
	v2 = ptr[1]/ptr[2] - v; v2 *= v2;
    } else {
	u2 = ptr[0] - u; u2 *= u2;
	v2 = ptr[1] - v; v2 *= v2;
    }

    dist = sqrt(u2 + v2);
    if (NEAR_ZERO(dist, 1.0e-4))
	return TRIM_ON;

    return TRIM_OUT;

}


/**
 * Process Case C curves;
 *
 * A check is placed here to determine if the u, v is on the curve.
 * Determine how many times the curve will cross the u, v axis. If the
 * curve crosses an odd number of times than the point is IN, else the
 * point is OUT. Since a Case C curve need processing a call to clip
 * the curve so that it becomes either Case B or Case A is required
 * to determine the number of crossings accurately. Thus we need to keep
 * the original curve and expect the calling routine to free the
 * storage. Additional curves are generated in this routine, each of
 * these new curves are processed, and then are deleted before exiting
 * this procedure.
 */
int
rt_process_casec(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
{

    struct edge_g_cnurb * clip;
    int jordan_hit;
    struct bu_list plist;
    int trim_flag = 0;
    int caset;

    /* determine if the u, v values are on the curve */

    if (rt_nurb_uv_dist(trim, u, v)  == TRIM_ON) return TRIM_IN;

    jordan_hit = 0;

    BU_LIST_INIT(&plist);

    if (nurb_crv_is_bezier(trim))
	rt_clip_cnurb(&plist, trim, u, v);
    else
	nurb_c_to_bezier(&plist, trim);

    while (BU_LIST_WHILE(clip, edge_g_cnurb, &plist)) {
	BU_LIST_DEQUEUE(&clip->l);

	caset = rt_trim_case(clip, u, v);

	trim_flag = 0;

	if (caset == CASE_B)
	    trim_flag = rt_process_caseb(clip, u, v);
	if (caset == CASE_C)
	    trim_flag = rt_process_casec(clip, u, v);

	nmg_nurb_free_cnurb(clip);

	if (trim_flag == TRIM_IN) jordan_hit++;
	if (trim_flag == TRIM_ON) break;
    }

    while (BU_LIST_WHILE(clip, edge_g_cnurb, &plist)) {
	BU_LIST_DEQUEUE(&clip->l);
	nmg_nurb_free_cnurb(clip);
    }

    if (trim_flag == TRIM_ON)
	return TRIM_ON;

    else if (jordan_hit & 01)
	return TRIM_IN;
    else
	return TRIM_OUT;
}


/**
 * This routine will be called several times, once for each portion of
 * the trimming curve. It returns whether a line extended from the
 * \<u, v\> point will cross the trimming curve an even or odd number of
 * times. Or the \<u, v\> point could be on the curve in which case
 * TRIM_ON will be returned. The algorithm uses the approach taken Tom
 * Sederburge and uses bezier clipping to produce caseA and caseB
 * curves. If the original trimming curve is a CASE C curve then
 * further processing is required.
 */
int
rt_uv_in_trim(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
{

    int quad_case;

    quad_case = rt_trim_case(trim, u, v);	/* determine quadrants */

    /* CASE A */
    if (quad_case == CASE_A)
	return TRIM_OUT;
    if (quad_case == CASE_B)			/* CASE B */
	return rt_process_caseb(trim, u, v);
    if (quad_case == CASE_C)			/* CASE C */
	return rt_process_casec(trim, u, v);

    bu_log("rt_uv_in_trim: rt_trim_case() returned illegal value %d\n", quad_case);
    return -1;
}


/**
 * This routines is used to determine how far a point is from the u, v
 * quadrant axes.
 *
 * Equations 3, 4, 5 in Sederberg '90 paper
 */
fastf_t
rt_trim_line_pt_dist(struct _interior_line *l, fastf_t *pt, int pt_type)
{
    fastf_t h;
    int h_flag;

    h_flag = RT_NURB_IS_PT_RATIONAL(pt_type);

    if (l->axis == 0) {
	if (h_flag) h = (pt[1] / pt[2] - l->o_dist) * pt[2]; /* pt[2] is weight */
	else h = pt[1] - l->o_dist;

    } else {
	if (h_flag) h = (pt[0] / pt[2] - l->o_dist) * pt[2];
	else h = pt[0] - l->o_dist;

    }

    return h;
}


/**
 * Return the SIGN of the value
 */
int
_SIGN(fastf_t f)
{
    if (f < 0.0)
	return -1;
    else
	return 1;

}


/**
 * We try to clip a curve so that it can be either Case A or Case C.
 * Sometimes one of the curves is still Case C though, but it is much
 * smaller than the original, and further clipping will either show that
 * it is on the curve or provide all Case B or Case A curves.  We try
 * to pick the best axis to clip against, but this may not always
 * work. One extra step that was included, that is not in the paper
 * for curves but is for surfaces, is the fact that sometimes the
 * curve is not clipped enough, if the maximum clip is less than .2
 * then we subdivide the curve in three equal parts, at .3 and .6 .
 * Subdivision is done using the Oslo Algorithm, rather than the other
 * methods which were prossed.
 */
void
rt_clip_cnurb(struct bu_list *plist, struct edge_g_cnurb *crv, fastf_t u, fastf_t v)
{
    fastf_t ds1, dt1;
    struct _interior_line s_line, t_line;
    int axis, i;
    fastf_t umin, umax;
    int coords;
    struct edge_g_cnurb * c1, *c2, *tmp;
    fastf_t m1, m2;
    int zero_changed;
    fastf_t *ptr;
    fastf_t dist[10];

    coords = RT_NURB_EXTRACT_COORDS(crv->pt_type);

    s_line.axis = 0;	s_line.o_dist = v;
    t_line.axis = 1;	t_line.o_dist = u;

    ds1 = 0.0;
    dt1 = 0.0;

    ptr = crv->ctl_points;


    /* determine what axis to clip against */

    for (i = 0; i < crv->c_size; i++, ptr += coords) {
	ds1 +=
	    fabs(rt_trim_line_pt_dist(&s_line, ptr, crv->pt_type));
	dt1 +=
	    fabs(rt_trim_line_pt_dist(&t_line, ptr, crv->pt_type));
    }

    if (ds1 >= dt1) axis = 0; else axis = 1;

    ptr = crv->ctl_points;

    for (i = 0; i < crv->c_size; i++) {
	if (axis == 1)
	    dist[i] = rt_trim_line_pt_dist(&t_line, ptr, crv->pt_type);
	else
	    dist[i] = rt_trim_line_pt_dist(&s_line, ptr, crv->pt_type);

	ptr += coords;
    }

    /* Find the convex hull of the distances and determine the minimum
     * and maximum distance to clip against. See the paper for details
     * about this step
     */

    umin = 10e40;
    umax = -10e40;
    zero_changed = 0;

    for (i = 0; i < crv->c_size; i++) {
	fastf_t d1, d2;
	fastf_t x0, x1, zero;

	if (i == (crv->c_size -1)) {
	    d1 = dist[i];
	    d2 = dist[0];
	    x0 = (fastf_t) i / (fastf_t) (crv->c_size - 1);
	    x1 = 0.0;
	} else {
	    d1 = dist[i];
	    d2 = dist[i+1];
	    x0 = (fastf_t) i / (fastf_t) (crv->c_size - 1);
	    x1 = (i+1.0) / (crv->c_size - 1);
	}

	if (_SIGN(d1) != _SIGN(d2)) {
	    zero = x0 - d1 * (x1 - x0)/ (d2-d1);
	    if (zero <= umin)
		umin = zero * .99;
	    if (zero >= umax)
		umax = zero * .99 + .01;
	    zero_changed = 1;
	}
    }

    if (!zero_changed)
	return;

    /* Clip is not large enough, split in thirds and try again */

    if (umax - umin < .2) {
	umin = .3; umax = .6;
    }

    /* Translate the 0.0-->1.09 clipping against the real knots */

    m1 = (crv->k.knots[0] * (1 - umin)) +
	crv->k.knots[crv->k.k_size -1] * umin;

    m2 = (crv->k.knots[0] * (1-umax)) +
	crv->k.knots[crv->k.k_size -1] * umax;

    /* subdivide the curve */
    c1 = (struct edge_g_cnurb *) nmg_nurb_c_xsplit(crv, m1);
    c2 = nmg_nurb_c_xsplit((struct edge_g_cnurb *) c1->l.forw, m2);

    tmp = (struct edge_g_cnurb *) c1->l.forw;
    BU_LIST_DEQUEUE(&tmp->l);
    nmg_nurb_free_cnurb(tmp);

    BU_LIST_INIT(plist);
    BU_LIST_INSERT(&c2->l, plist);
    BU_LIST_APPEND(plist, &c1->l);
}


int
nmg_uv_in_lu(const fastf_t u, const fastf_t v, const struct loopuse *lu)
{
    struct edgeuse *eu;
    int crossings=0;

    NMG_CK_LOOPUSE(lu);

    if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
	return 0;

    for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
	struct edge_g_cnurb *eg;

	if (!eu->g.magic_p) {
	    bu_log("nmg_uv_in_lu: eu (%p) has no geometry!!!\n", (void *)eu);
	    continue;
	}

	if (*eu->g.magic_p != NMG_EDGE_G_CNURB_MAGIC) {
	    bu_log("nmg_uv_in_lu: Called with lu (%p) containing eu (%p) that is not CNURB!!!!\n",
		   (void *)lu, (void *)eu);
	    bu_bomb("nmg_uv_in_lu: Called with lu containing eu that is not CNURB!!!\n");
	}

	eg = eu->g.cnurb_p;

	if (eg->order <= 0) {
	    struct vertexuse *vu1, *vu2;
	    struct vertexuse_a_cnurb *vua1, *vua2;
	    point_t uv1, uv2;
	    fastf_t slope, intersept;
	    fastf_t u_on_curve;

	    vu1 = eu->vu_p;
	    vu2 = eu->eumate_p->vu_p;

	    if (!vu1->a.magic_p || !vu2->a.magic_p) {
		bu_log("nmg_uv_in_lu: Called with lu (%p) containing vu with no attribute!!!!\n",
		       (void *)lu);
		bu_bomb("nmg_uv_in_lu: Called with lu containing vu with no attribute!!!\n");
	    }

	    if (*vu1->a.magic_p != NMG_VERTEXUSE_A_CNURB_MAGIC
		|| *vu2->a.magic_p != NMG_VERTEXUSE_A_CNURB_MAGIC)
	    {
		bu_log("nmg_uv_in_lu: Called with lu (%p) containing vu that is not CNURB!!!!\n",
		       (void *)lu);
		bu_bomb("nmg_uv_in_lu: Called with lu containing vu that is not CNURB!!!\n");
	    }

	    vua1 = vu1->a.cnurb_p;
	    vua2 = vu2->a.cnurb_p;

	    VMOVE(uv1, vua1->param);
	    VMOVE(uv2, vua2->param);

	    if (RT_NURB_IS_PT_RATIONAL(eg->pt_type)) {
		uv1[0] /= uv1[2];
		uv1[1] /= uv1[2];
		uv2[0] /= uv2[2];
		uv2[1] /= uv2[2];
	    }

	    if (uv1[1] < v && uv2[1] < v)
		continue;
	    if (uv1[1] > v && uv2[1] > v)
		continue;
	    if (uv1[0] <= u && uv2[0] <= u)
		continue;
	    if (ZERO(uv1[0] - uv2[0])) {
		if ((uv1[1] <= v && uv2[1] >= v) ||
		    (uv2[1] <= v && uv1[1] >= v))
		    crossings++;

		continue;
	    }

	    /* need to calculate intersection */
	    slope = (uv1[1] - uv2[1])/(uv1[0] - uv2[0]);
	    intersept = uv1[1] - slope * uv1[0];
	    u_on_curve = (v - intersept)/slope;
	    if (u_on_curve > u)
		crossings++;
	} else
	    crossings += rt_uv_in_trim(eg, u, v);
    }

    if (crossings & 01)
	return 1;
    else
	return 0;
}


/*
 * Local Variables:
 * mode: C
 * tab-width: 8
 * indent-tabs-mode: t
 * c-file-style: "stroustrup"
 * End:
 * ex: shiftwidth=4 tabstop=8
 */
